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Abstract 



Networks of coupled neural systems represent an important class of models in com- 
putational neuroscience. In some applications it is required that equilibrium points 
in these networks remain stable under parameter variations. Here we present a gen- 
eral methodology to yield explicit constraints on the coupling strengths to ensure 
the stability of the equilibrium point. Two models of coupled excitatory-inhibitory 
oscillators are used to illustrate the approach. 
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1 Introduction 



We consider neural networks of the form 



N 

x* = F(x*) + £G,,H(x^ (1) 
j'=i 

where x* is the M-dimensional state vector of the ith node. Each node can 
either be a single neuron (M = 1) as in Hopfield types of models [13,27], or a 
group of neurons (M > 1), representing e.g. the cortical column of interacting 
excitatory and inhibitory neurons [20,26,29]. The dynamics of the individual 
node is given by x* = F(x l ) and H : R M — > R M is the coupling function. The 
coupling matrix is G = [Gy] where gives the coupling strength from node 
j to node %. 

Without loss of generality assume that the origin is a stable equilibrium point 
for the individual node and remains an equilibrium point for the network. The 
stability of the origin under coupling strength variations is the main concern 
of the present work. This problem is mainly motivated by some computational 
considerations. For example, a class of models assert that the background state 
of the network, represented by the equilibrium point at the origin, should 
be quiescent in the absence of input [2,3,4,6,17,18,19,28,31]. External inputs, 
treated as a slowly increasing and then decreasing function of time, can lead 
the network through a Hopf bifurcation to an oscillatory state and then return 
it to its background or equilibrium state once the input has been removed. 
This natural reset mechanism, requiring the origin to be a stable equilibrium 
point, makes the network ready for the next computational cycle. To endow 
the oscillatory network the ability to differentiate patterns of inputs, statistical 
learning takes place wherein the coupling strengths between the network units 
change according to certain learning rules. Without careful consideration the 
learning related parameter changes can potentially alter the stability of the 
background state, thereby defeating the computational picture established 
earlier. It is thus desirable to have constraints on the individual coupling 
strengths that can be incorporated into the learning rules so that the stability 
of the equilibrium point is ensured for all time. 

Previous work on stability constraints have mainly concentrated on recurrent 
networks of the Hopfield type [1,5,7,8,11,13,14,16,22,23,24,25,27,30] with M = 
1. In this paper we consider a general approach that leads to stability bounds 
on the individual coupling strengths in recurrent networks with more complex 
local dynamics. Two explicit models of coupled neural populations will be 
used to illustrate our approach. 
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2 Theory 



Our approach consists of three steps. 

Step 1. For simplicity, let F(0) = 0, H(0) = 0, and the real parts of the 
eigenvalues of the Jacobian DF(0) be negative so that the origin is stable for 
the individual node. 

Linearizing Eq. (1) around the origin gives (in matrix form) 



where S = (x^x 2 , • • • ,x JV ). According to the Jordan canonical form theory, 
the stability of Eq. (2) is determined by the eigenvalue A of G. Let the cor- 
responding eigenvector from G T be e and let u = Se. The equation for u 
reads 



The origin of Eq. (1) is stable if this equation is stable for all the eigenvalues 
of G. This is true even when the coupling matrix is defective [12]. 

Step 2. To proceed further we treat A in Eq. (3) as a complex control parame- 
ter. Denote by Q the region in the Re(A)- Im(A) plane where all the eigenvalues 
of (DF + A • DH) have negative real parts. Clearly, the equilibrium point is 
stable if all eigenvalues of G lie within Q. We henceforth refer to Q as the sta- 
bility zone. A schematic of Q is shown in Figure 1. We note that Q is usually 
obtained numerically. For some situations analytical results are possible (see 
below). 

Step 3. Thus far the stability criteria are stated in terms of the eigenvalue 
of G. The goal in this work is to directly constraint the coupling strengths 
themselves. This is done by making use of the Gershgorin disc theorem [15]. 

Given an n x n matrix A = [a^], the Gershgorin theorem states that all 
eigenvalues of A are located in the union of n discs (called the Gershgorin 
discs) where each disc is given by 

{z e C : \z - au\ < \ a ji\}, i = l,2,...n. 



S = DF-S + DH-S-G T , 



(2) 



u = [DF + A ■ DH]u. 



(3) 



Alternative forms of the n discs are [15]: 
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{z e C : \z - au\ < \ a ij\}, i = l,2,...n. 
Combining the two, we have the form used in the remainder of this paper: 

{z e C: \z-aa\ < ^^2(\aji\ + \aij\)}, i = l,2,...ra. (4) 

This form is more intuitive since it involves incoming and outgoing coupling 
strengths for a given node. 

The stability conditions for the equilibrium point can now be stated as follows: 

(1) The center Gu (i — 1, 2, . . . , N) of every Gershgorin disc of G lies inside 
the stability zone Q; 

(2) The radius of every Gershgorin disc is shorter than the distance from the 
center of the disc to the boundary of fi. 

In other words, letting S(x) denote the distance from point x on the real axis 
to the boundary of Q, stability of the equilibrium point is ensured if 

(G ri ,0)en and ^£(\ G ji\ + \Gij\) < S(G U ) (5) 
for % = 1,2,..., N. 



3 Examples 

3.1 The case of M = 1 

When one dimensional systems are coupled together, the matrices DF and DH 
are reduced to real numbers. Representing them by fi and v respectively, the 
stability zone is easily obtained as Re(A) < —^jv. The distance from the center 
of the ith Gershgorin disc to the boundary of Q is given by 5(Ga) = —fi/u—Ga. 
Using Eq. (5) we obtain the stability conditions as 

^E(l^l + l^'l) + ^<-^- (6) 

This result was obtained before in [13,27]. 
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3.2 A coupled oscillator model with M = 2 



The general topology for the model is shown in Figure 2. The basic unit in the 
model is a neural population consisting of either excitatory or inhibitory cells 
[2,17,21,29]. The functional unit in the network is a cortical column consisting 
of mutually coupled excitatory and inhibitory populations. The columns are 
then coupled through mutually excitatory interactions to form the network. 

A single column is described by two first order differential equations 



dx 

— + ax = -k ei Q(y, Q m ) + I, 
^jt + by = k ie Q(x,Q m ). 



(7) 



Here x, y represent the local field potentials of the excitatory and inhibitory 
populations, respectively, and / is the input (I — in the subsequent analy- 
sis). The constants a, b > are the damping constants. The parameter k ie > 
gives the coupling gain from the excitatory (x) to the inhibitory (y) population 
whereas k e i > represents the strength of the reciprocal coupling. The non- 
linear neuronal interaction is realized through the sigmoid function Q(-,Q m ) 
where Q m is a parameter controlling the slope of the function. Here we only 
need to specify that Q(0, Q m ) = and Q'(0, Q m ) = 1. 

The iV columns are coupled together in the following fashion: 



dx n 
~dt 

dy n 
dt 



N 



+ ax n k e iQ(jj n , Qm) + y ' c n pQ(^Xp, Qm) -\- I n , 



P =i 



(8) 



"I - by n ki e Q{x n) Qm)i 



where the columns are indexed byn = 1,2, ... ,N and the coupling strength 
c np is the gain from the excitatory population of column p to the excitatory 
population of column n. 

Variables used in Eq. (3) can be explicitly evaluated for the present model as 



where we have used the fact Q'(Q, Q m ) = 1. 



6 



To discover the stability zone we study the eigenvalue a of the matrix (DF + 
A • DH) as a function of A. The characteristic polynomial of this matrix is 
given by 

f(a) = a 2 + a(a + b — A) + (k ei k ie + ab — bX). 
For an arbitrary coupling matrix G, its eigenvalues A could be complex: 

X = X R + iXj. 
Then the characteristic polynomial becomes 

f(a) = a 2 + a(a + b — X R — iXi) + (k ei k ie + ab — bX R — ibXi). 

The range of parameter values which gives Re(a) < can be determined by 
applying the generalized Routh-Hurwitz criterion (see Appendix I). Following 
this procedure, consider —if(ia): 

— if(ia) = ia 2 + a(a + b — Xr) — iaXi — i(k ei k ie + ab — bXn) — bXj. 
This has to be put into the following standard form: 

— if(ia) = b a 2 + bict + b 2 + i[a a 2 + a\(x + a 2 ). 
Comparing the two equations we get 

a = l, a 1 = -X I , a 2 = -(k ei k ie + ab - bX R ), 
b = 0, bi = (a + b - X R ), b 2 = -6A/. 

Applying the generalized Routh-Hurwitz criterion, we have Re (a) < if the 
following two conditions are met: 



1 -A, 
(a + b - Xr) 



> 



and 
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1 -A/ -(k ei k ie + ab- bX R ) 

0{a + b-X R ) -6A/ 

1 -A/ -{keihe + ab - b\ R ) 

(a + fe-Aij) -6A/ 



> 0. 



Evaluating the above determinants and simplifying, we get 



(a + b- X R ) >0, 

(/c ei /c je + ab - bX R )(a + 6 - A^) 2 - bX 2 j(X R - a) > 0. 



(9) 



Solving the inequalities, the stability zone Q (see Figure 3) is found to be the 
region to the left of the curve 



, _ (keihe + ab- bX R )(a + b- X R ) 2 



b(X R - a) 



(10) 



The pointed tip of the curve in Figure 3 along the real axis is given by (min(a+ 
6, a + k ie k e i/b), 0) and it corresponds to the symmetric coupling case. 

The distance 5(Ga) from the center of the ith Gershgorin disc to the boundary 
is (see Appendix II for more details) 



5 (Go) = y (a - Gu) 2 - b 2 - 2k ie kei + 2yf k ie k ei [2b(a + b- Go) + k ie k e 
So the stability conditions [cf. Eq.(5)] are given by 



Gu < min(a + b, a + k ie k ei /b), 



a - Gu) 2 - b 2 - 2k ie kei + 2y k ie k e i[2b(a + b- Gu) + k ie k e 



We note that, since the boundary curve of the stability zone asymptotically 
approaches the straight line X R = a, we can use this line to define a new 
stability zone to obtain some simpler stability constraints. The distance to 
the new boundary is easily found to be 



Si — | a — Gu I . 
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In this case, the stability condition simplifies to 



^EO^I + fail) +Gu<a, i = 1, 2, . . . , N. (11) 

This simplified condition is a good approximation if min(a + b, a + ki e k e i/b) is 
sufficiently close to a. We further note that Eq. (11) is satisfied if 

\Gij\<a/N, i,j = 1,2,..., AT. 
That is, the equilibrium point is stable if 



\c np \ < a, V n,p = 1,2, . . ., N. 

This simple stability bound on the individual coupling strengths can be very 
useful in practice. 



3.3 A coupled oscillator model with M = 4 



The previous model represents a neural population by a first order differential 
equation. This has the property that its impulse response has a instantaneous 
rise phase. Here we consider another model where the neural population is a 
second order differential equation possessing a finite rise and decay impulse 
response. Each individual column is described by a system of two second order 
differential equations [9]: 



d^x dx 

— + (a + b)— + abx = -k ei Q{y, Q m ) + I, 
Cpy cly 

— + (a + b)— + aby = k ie Q{x, Q m ). 



(12) 



The parameters have the same interpretation as before. The N column equa- 
tions are given by: 
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+ (a + 6)^jT + afa ™ = -keiQiVn, Qm) + 



A' 



A r 



^ ' CnpQ{%p) Qm) I-ni 



(13) 



P =i 



+ (a + b)—^ + aby n = k ie Q(x n , Q m ), 
where the same network topology in Figure 2 applies. 

We first consider the stability of the single column equations given in Eq. 
(12). When the input / is zero, the origin x = 0,y = is an equilibrium 
point. In order to study its stability properties, we convert the above second 
order differential equations to the following system of first order differential 
equations: 



— = -(a + b)z 2 - abzi - k ei Q(z 3 , Q m ), 



dzi 
~dt 

dz 2 

~dl 

dz 3 

dt 

dz4 

~dt 



Z2, 



z 4 , 



— = -(a + b)zi - abz 3 + k ie Q(z 1 , Q m ), 



where 



dx 



dy 



Z\ = X, Zo = — , Zs — V, Z4 = —. 

dt' y dt 



The Jacobian matrix DF is obtained as 



DF 



( 
-ab 





1 

(a + b) 







— k ■ 


—ab 





1 

-(a + b)J 



(14) 



Here we have used the fact that Q'(0, Q m ) = 1- For stability of the origin, the 
real parts of all eigenvalues of DF should be less than zero. The eigenvalues 
are determined from the characteristic equation: 



A + 2(a + b)X 6 + (a 2 + Aab + b 2 )\ 2 + 2(a 2 b + ab 2 )\ + k ie k ei + a 2 b 2 = 0. 

Applying the Lienard-Chipart criterion (see Appendix I), the real parts of all 
eigenvalues are negative if the following inequalities be satisfied: 
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a 2 b 2 + k ie k ei > 0, 
a 2 b + ab 2 > 0, 
a + b > 0, 
—k ie k ei + ab(a + b) 2 > 0. 

Since a, 6, k e i, k{ e > 0, the first three inequalities are automatically satisfied. 
After simplification, the last inequality can be written as: 

kiek e i < ab(a + b) 2 . (15) 

To summarize, the origin is stable for the single column equations if the above 
condition is satisfied. Henceforth, we will assume that this is true. 

Next, we consider the stability of a network of coupled columns given in Eq. 
(13). Here 

\S*\ny — j 

and 



DH = 



(0 








o\ 


1 























Vo 








o) 



As before, we examine the eigenvalue a of the matrix DF+X-DH as a function 
of A. The characteristic polynomial of this matrix is given by 

f(a) = a 4 + 2(a + b)a 3 + [(a + b) 2 + 2ab - \]a 2 

+ [2ab(a + b) - X(a + b)]a + [a 2 b 2 - abX + k ie k ei ]. 

For complex A, we are not able to obtain an analytical form for the stability 
zone fl, since the characteristic equation results in a 8th order polynomial 
when applying the generalized Routh-Hurwitz criterion. However, numerical 
results are always possible. Figure 4 shows the stability zone Q when a = 
0.22, b = 0.72, k ie = 0.1, k ei = 0.4. After numerically finding the distance 
5 (Go) from the center of the ith Gershgorin disc to the boundary curve, Eq. 
(5) can again be used to give the stability criteria. 

If the coupling is symmetric, which implies that A is real, the stability bound- 
ary is just the rightmost tip of the curve along the real axis in Figure 4. Then 
the distance 5 is given by the absolute difference between the coordinates 
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of the tip point and the center of the ith Gershgorin disc. This tip can be 
determined as follows. 

Again applying the Lienard-Chipart criterion (see Appendix I), the real parts 
of all eigenvalues are negative if the following inequalities are satisfied: 



a 2 b 2 — abX + k ie k e i > 0, 
2ab(a + b) - A(a + 6) >0, 
(a + 6) >0, 

A 2 - 2(a + b) 2 \ + 4(a 3 6 + 2a 2 b 2 + ab 3 - k ie k ei ) > 0. 



(16) 



Since a, b are positive, the third inequality is automatically satisfied. After 
simplification, the first two inequalities become: 



ki e k e i + a 2 b 2 



ab 



A < lab. 



The last inequality is of the form 

aiA 2 — a 2 \ + a 3 > 0, 
where 



ai = 1, a,2 = 2(a + 6) 2 , as = A[ab(a + b) 2 — k ie k e i\. 

Note that ai, a 2 are obviously positive. It turns out a 3 is also positive because 
of the local stability condition derived in Eq. (15). The quadratic function 
aiA 2 — a 2 A + a 3 with a±, a 2 , a 3 positive has a unique global minimum at 
A = a 2 /2a 1 . Thus the minimum occurs at a positive value of A. It is also seen 
that 



a\ - 4aia 3 = 4 (a + 6) 4 - A[ab(a + b) 2 - k ie kei] 
This can be simplified as 



4aia 3 = 4 (a — b ) + Ak ie k, 



which is positive since k ie k e i is positive. Thus both the zeros of the quadratic 
function (we will denote them r\\ and r/ 2 with r\\ < r] 2 ) are real. Further, since 
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a3 > and the global minimum occurs at a positive value, 7/2 > rji > 0. 
Consequently, the last inequality is satisfied when A < 771 or A > r] 2 where 



771,2 = (a + b) 2 ± V(a 2 - 6 2 ) 2 + Ak ie k ei . 

Note that 771 is explicitly seen to be positive by applying Eq. (15). Further, 
772 > (a + b) 2 > lab. Thus the inequality A > 7/2 > 2a6 is not possible given 
the stability condition A < 2ab derived earlier. Therefore the last inequality 
in Eq. (16) reduces to A < 771. 

Summarizing, we get the following set of stability conditions: 

ki e k e i + a 2 b 2 

A < 1 ' 

ab 

A < 2ab, 
A<?7i. 

fc. . _|_ Q 2 ^ 2 

Let re = mini %e e% , 2ab, rjA, then all these inequalities will be simul- 

ab 

taneously satisfied if 

A < re. (17) 



Thus the rightmost tip of the boundary curve along the real axis is (re, 0). 
Therefore the distance function 5(Gu) is given by 



5(G U ) = \k-G u \, 1 = 1,2,..., TV. (18) 

Applying Eq. (5), we obtain the following stability condition for the present 
model with symmetric couplings: 



^Ed^l + \ Q ij\) + Gu<k } % = 1, 2, . . . , N. (19) 

As we discussed before, this condition is satisfied if the individual coupling 
strengths obey the following stability constraints: 



I dp I < re, for c np = c pn , n,p = 1,2, . . . , N. (20) 
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4 Conclusions 



We have presented a general method for studying the stability of the equi- 
librium state in neural network models. When the single-neuron coupled net- 
works, such as Hopfield type of models, are studied, the stability result from 
our general approach coincides with the known result found in the literature. 
As a harder application, two typical neural population models where the in- 
dividual nodes are higher dimensional were considered. The stability of the 
first model, a coupled network of two dimensional systems, was solved com- 
pletely. For the second model, a coupled network of four dimensional systems, 
stability criteria for symmetric coupling was given analytically. For the non 
symmetric case, our method was used to obtain numerical criteria. Through 
the above examples we have demonstrated that our general method is appli- 
cable to arbitrary neural networks where the individual nodes can themselves 
be high dimensional. When the dimension of the individual node is not too 
high, analytical results are possible. 

iFrom the stability criteria, we also derived simple bounds on the coupling 
strengths which ensure stability. These bounds put a limit on the magnitude 
of change that the coupling strengths can undergo in the process of statistical 
learning. 
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Appendix I 



In this Appendix, we state the Lienard-Chipart and generalized Routh-Hurwitz 
criteria. The statements are taken directly from Gantmacher [10] and are given 
here for the sake of completeness. 



A. Lienard-Chipart Criterion 



Consider a real polynomial 



f(z) = a z n + a x z n - x + 



+ On, 



with a > 0. Necessary and sufficient conditions for all the zeros of the poly- 
nomial to have negative real parts can be given in any one of the following 
forms [10]: 

(1) a n > 0, a„_ 2 > 0, . . .; Ai > 0, A 3 > 0, . . ., 

(2) a n > 0, a„_ 2 > 0, . . .; A 2 > 0, A 4 > 0, . . , 

(3) a n > 0, a„_i > 0, a„_ 3 > 0, . . .; A 1 > 0, A 3 > 0, . . ., 

(4) a n > 0, a„_i > 0, a„_ 3 > 0, . . .; A 2 > 0, A 4 > 0, . . .. 

Here A p is the Hurwitz determinant of order p given by the formula 



0\ fl 3 CI5 . . . 

CLq 0,2 CL4 ... 

ai a 3 . . . 

do «2 «4 



, p=l,2,...,n, 



where Ok = for k > n. In the literature, the equivalent Routh-Hurwitz 
criterion is usually used. But the Lienard-Chipart is better since the number 
of determinants that have to be evaluated is half the number that have to 
be evaluated for the Routh-Hurwitz criterion. This leads to a simpler set of 
inequalities that need to be evaluated. In the main text, we use the third form 
of the Lienard-Chipart criterion given above. 
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B. Generalized Routh-Hurwitz Criterion 



Consider a polynomial f(z) with complex coefficients. Suppose that 



f(iz) = b z n + b x z n - x + ■ ■ ■ b n + t(a z n + a lZ n ~ l + ■ ■ ■ + a n ), 

where a , a±, . . ., a n , b , bi, . . ., b n are real numbers. If the degree of f(z) is 
n, then b + ia Q ^ 0. Without loss of generality, we may assume that a Q ^ 0. 
Otherwise, we consider the polynomial g(z) = —if(z) and repeat the analysis 
for this polynomial. Both f(z) and g(z) have the same set of zeros and so no 
information is lost. This is the case considered in the main text. 

If V2« 7^ 0, then all the zeros of f(z) have negative real parts if 



V 2 >0, V 4 >0, V 2 n>0, 
where 



2p 



a ai 

bo h 

ao 

b 



a 2p-l 
blp-2 



p= 1,2,..., n, 



where = bk = for k > n. Note that the condition V2n 7^ would be 
satisfied for a generic set of parameter values. This is especially true in our 
case where a^, bk are functions of system parameters. 



Appendix II 



The distance 7 from the center (Gu, 0) of the ith Gershgorin disc to any point 
on the boundary of the stability zone is given by 
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(A_r — Go) 2 + \] 



Substituting A/ from Eq. (10) and differentiating with respect to A/, we have 
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drf_ _ 
dX R 

(a + b-X R ) 2 [(Xr- a) 2 - b 2 ](ab + k ie k ei - bX R ) 

— LrjjJ 77 v 1 777 To ' 

(X R - a) b(X R - a) 2 



dl 2 

Setting — — = 0, we get two solutions: 
dX R 



X R = a±bi 



k- h ■ 



2b(a + b-Gu + k ie k ei ) ' 



Since the boundary of Q lies to the right of the point (a,0), we can discard 
the smaller solution. Substituting the remaining solution in the equation for 
7 2 and taking the square root, we get the shortest distance as: 



H Imin 

\J(a- Go) 2 - b 2 - 2k ie k ei + 2\Jk ie k ei [2b(a + b - Go) + k ie k ei ], 
i = l,2,...,N. 
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Figure Caption 

Figure 1: Schematic of the stability zone. 
Figure 2: Schematic of the network configuration. 
Figure 3: Stability zone for model Eq.(8) 
Figure 4: Stability zone for model Eq.(13) 
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